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ABSTRACT 



Application of classical windows to time series data is a means of enhancing the 
performance of the periodogram. Use of these classical windows results in the broadening 
of the spectral mainlobe. A Kalman filter will smooth spectral data by dividing the 
periodogram of unwindowed time series data into piecewise constant segments, ideally into 
noise-only and signal-only segments. This allows for a more accurate representation of the 
mainlobe of the original periodogram. The Kalman filter was modified to alter the filter 
parameter /3 during filtering to provide maximum smoothing during the noise-only segment 
and maximum sensitivity in the vicinity of the spectral peak of the periodogram. This 
modification results in a smoothing of the noise-only portion of the periodogram while 
leaving the main spectral peak essentially unaltered. A second modification was made to 
substitute the original raw values of the periodogram for the filter estimates during detected 
up-transitions while smoothing the noise-only segments in the same manner as in the original 
Kalman filter algorithm. This further enhances the preservation of the mainlobe of the 
periodogram and lowers the noise floor 1 to 3 dB over that of the original Kalman filter. 
These processes were further enhanced by stacking the output periodograms and displaying 
them as LOFAR output on the Sun workstation. NCAR graphics grey-toning subroutine is 
used to generate the LOFAR displays. 
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I . INTRODUCTION 
A. MAN-MACHINE INTERFACE 

As the computational powers of the computer continues to 
improve, the ability to process large amounts of data in 
relatively short periods of time has also improved. In 
addition it is now possible do much of these computations on 
desk-top computers in the office or even at home. In the past 
the results (i.e., output) of these computations were provided 
in the form of paper printout. The individual who must digest 
this information and make a decision or draw conclusions from 
this data was faced with a serious problem. The value of the 
information provided will depend on the ability of the 
decision maker to understand what is being presented. 

Recent advances have made the computer very efficient at 
gathering and storing data. The problem is retrieving the 
data in a usable form and presenting it to the user (i.e., 
decision maker) . [Ref. 1] 

The use of computer graphical displays is an effective 
solution to this problem. Human information processing 
operations depend on the sensory reception of relative 
stimulation. The information we are interested in is 
contained in this stimuli. This stimuli may be coded in 
visual displays and reproduced. The reproductions can be 
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modified by enlargement, amplification, filtering or 
enhancement. When indirect sensing a; -plies, there is a human 
factors aspect that must enter into the design of these 
displays. [Ref. 2]. 

The human brain has a vast capacity for non-verbal 
perception. This trait can be optimized since computers have 
the ability to graphically communicate up to 100,000 times 
more effectively than statistical printouts. It is estimated 
that the total capacity of the visual channel can be as high 
as 30-40 million bits of information per second. This is 
equivalent to 48-72 million words per minute, assuming seven 
bits per character and five characters per word. Since humans 
can not read this fast the effectiveness of graphics over 
alpha-numeric printouts is clear [Ref. l:p.5]. 

Often people can sense many sources of information 
directly without much problem. There are however, many 
instances where direct sensing is not adequate. Some of these 
situations are mentioned below: 

1. Information is at or below threshold values that may need 
to be amplified by some means (i.e., electronic, 
visually, etc...). 

2 . Information may be to large or massive and require 
reduction or consolidation. 

3. Information could be embedded in noise that may need to 
be filtered out. 

4. Stimuli may be of the type that could be better sensed if 
converted to some other type of sensory display (e.g., 
graphs or color coding to represent quantitative data) . 
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Generally speaking, the information being presented by 
displays is either dynamic or static. Dynamic information is 
continually changing or may change with time. The randomness 
associated with the noise in which a signal may be present is 
considered dynamic. Displays may also present static 
information. This could be information that does not change 
with time of for information that remains constant for some 
fixed period of time. We are interested in a dynamic display 
in this discussion. For a detailed listing of specific 
examples of dynamic and static displays see [Ref. 2:p.49]. 

Selecting the proper display is very important and will 
depend on the type of data being evaluated. Some general 
guidelines for the selection of visual displays follows: 

1. Use the simplest display concept that will most rapidly 
provide the user with the information needed to properly 
interpret the data of interest. The more sophisticated 
and complicated the display the longer it will take the 
user to digest and the more likely the user will be to 
misinterpret the data being displayed. 

2. Use the least precise format needed to covey the desired 
information. Overemphasizing the accuracy of readout 
required can result in increasing the users response 
time, add to users fatigue or mental stress, and can 
ultimately cause unnecessary mistakes. 

3. Use a display that is natural to the user. The picture 
the user sees should convey the expected interpretation 
that is associated with the data being displayed. 

4. Use the most practical display technique for the expected 
viewing environment and operating viewing conditions. 
[Ref. 3]. 

Consider the case where the data of interest is a "n by n" 
matrix of considerable size and the user is interested in 
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locating the maximum value in that matrix while also 
recognizing any consistencies between the maximum values in 
each row. In addition, suppose that the exact numerical value 
of those maximum numbers are less important than their 
locations in their respective rows. It would be useful if the 
maximum value of the matrix could be assigned a visual 
parameter and all the other numerical values could be assigned 
a visual parameter that could be scaled relatively to that 
maximum. This visual coding takes up significantly less space 
than would be required when displaying the numerical value of 
each piece of data and would allow one look at the entire 
matrix at once. Consider a coding system in which the visual 
scale uses grey tones over the range of white to black. The 
largest value in the data set is coded to represent the 
darkest shade and the smallest value is assigned the absence 
of color (i.e. white) . All values in between will be assigned 
grey-tones relative to their numerical standing in 
relationship to the maximum and minimum values. This is the 
basic idea behind the grey-tone (LOFAR) display utilized in 
this research. 

B . RESEARCH OVERVIEW 

The estimate of the spectral content of a segment of a 
process is most commonly accomplished by the use of a 
periodogram. The performance of the periodogram has 
significant shortcomings when dealing with the detection of 
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signals in noise. One of these shortcomings is the window 
function sidelobe effects that result from processing data 
sets of finite length. By using tapered window functions, 
such as the Hamming, the von Hann, or the Bartlett window, one 
can improve the performance of the periodogram by minimizing 
the effects caused by the discontinuities of the data at the 
boundaries [Ref. 4]. The periodogram also has a large 
variance (independent of data length used). Typically, the 
mean equals the standard deviation in the noise-only case. An 
other method of improving the performance of the periodogram 
is by averaging over a succession of sequential periodograms . 
The smoothing of the spectral estimate is due to the reduction 
of the variance of the estimate. Loss of resolution and 
broadening of the main spectral peaks almost always results 
from the use of either of these techniques. In the first 
case, the broadening is caused by the convolution of the 
window Fourier transform with the Fourier transform of the 
data, while in the second case the broadening potentially 
comes from having to shorten transform lengths. William Go 
[Ref. 5] performed research that examined the 
application of a Kalman filter to a rectangular windowed data 
set. In his work he was able to demonstrate that the Kalman 
filter could filter (i.e. smooth) the noise portions of the 
spectral estimate while leaving the main spectral peaks 
essentially unaltered. The results realized were preservation 
of the original periodograms resolution while accentuating the 
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dominant spectral peaks as they rise out of the noise floor. 
Go used single and multiple sinusoids in Gaussian white noise 
to evaluate the Kalman filter's performance in signal 
detection and resolution at different input signal-to-noise 
ratios (SNR) for multiple noise realizations. He also 
examined the effects of altering the filter's detection 
parameter (called beta, jS) and the data /transform lengths. It 
was demonstrated by Go that the Kalman filter did appear to 
out perform the Hamming window in signal detectability as well 
as to provide better preservation of the spectral peaks. It 
is our purpose to take a more in depth look at the Kalman 
filter used in Reference 4. 

1. Filter Pareuneter (/3) 

The filter parameter is currently chosen by the user 
in an experimental fashion. The larger the /3 the less 
sensitive the filter is to the dynamics of the spectrum and 
the smaller the /S the more sensitive the filter will behave. 
There are some obvious trade-offs here. In [Ref. 5], a 
/ 3=3 00000 appeared to be the most effective for detecting the 
spectral peak of signals with FFT output SNR's of 9 dB or 
greater. The purpose here is to attempt to refine the choose 
of for signal detection of signals with FFT output SNR's 
less than 9 dB. 
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a . Dynamic Beta 

An area of interest is the study of the effects of 
altering /S during the filtering process. If the algorithm is 
altered so that once a detection occurs the value of /3 is 
reduced to make the filter very sensitive (thus preserving the 
original shape of the periodograms) during the down transition 
of the spectral peak, then it may be possible to retain the 
resolvability of two closely spaced spectral components. 
Similarly, once past a given spectral component, /3 is reset to 
its original value and the filtering continues. 

b. Modified Kalman Filter 

Since we know that spectral components will appear 
as positive peaks, and usually are of relatively short 
duration, the Kalman filter algorithm was modified to replace 
the filter estimate with the original periodogram data while 
the filter is detecting up transitions in the data. 

2 . Output Presentation 

For very large signal lengths it is cumbersome to look 



at the periodograms line 


by 


line 


(i.e. transform 


by 


transform) . 


It 


would be 


useful to 


somehow stack 


the 


periodograms 


and 


look at 


the 


final 


outcome of all 


the 



transforms in their entirety such as in a 3-Dimensional plot 
or intensity plot (LOFAR) . An additional objective of this 
research is the development of a display of these stacks of 
periodograms in a grey-tone image on the Sun workstations. 
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The graphics package used in this research is NCAR graphics. 
The LOFAR outputs are created using the subroutine HALFTON 
which draws a halftone picture from data stored in a 
rectangular array. These arrays are read in row by row. The 
subroutine assigns to the largest value in the data the 
darkest grey-tone (i.e. highest intensity), and conversely 
assigns the lightest grey-tone to the lowest value in the 
data. Sixteen intensity levels for the mapping of data can be 
used in a halftone image. To each intensity corresponds an 
equal range in the data (i.e. linear mapping) . Appendix D is 
a step-by-step users guide that provides the necessary 
guidance to produce LOFAR output on a NFS Sun work station. 
The NFS Sun stations use the UNIX operating system via an 
Internet Networking System. 
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II 



SPECTRAL ESTIMATION 



A. BACKGROUND 

Fourier analysis has been an analytical tool for use with 
analog signal processing for continuous-time signals. Due to 
the vast improvements made in the computational power of the 
digital computers, the old school approach is no longer the 
only method of choice. There are many applications today in 
analog signal processing where it is preferred to sample the 
original signal, process it on a digital computer utilizing a 
discrete-time system, then convert it back to an analog signal 
if desired. One technique to accomplish this is the discrete 
Fourier analysis, which is the discrete version of continuous- 
time Fourier analysis. Fourier-based estimation, also called 
classical spectral estimation, includes two main nonparametric 
methods. The first is the periodogram method, which has 
limited resolution and a variance of the spectral estimate 
that remains constant independent of increases in data length. 
Windowing techniques may be used to improve the variance. 
These techniques are extremely efficient and produce adequate 
results for many dissimilar types of signals. The second 
Fourier-based method is the Blackman-Tukey method which 
utilizes the Fourier transform of an estimate of the 
autocorrelation function. A problem of the Fourier-based 
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methods is caused by the windowing of the data during 
processing. Although windowing minimizes the effects caused 
by discontinuities of data at the boundaries, a drawback with 
the use of windows is that leakage in the spectral domain 
nearly always results. This is where energy in the main lobe 
of the spectral response leaks into the nearby sidelobes, 
which can obscure and misrepresent other spectral frequency 
components that may be present. Using a non-tapered window 
(i.e., rectangular window) provides the best resolution of any 
window currently in use although the spectral leakage 
resulting from using a rectangular window is worse than 
realized with any other window. [Ref. 6]. 

B . PERIODOGRAM 

There are many procedures used to estimate the spectral 
content of a sample of a process. The popular choice remains 
the periodogram primarily due to its ability to economically 
and effectively be implemented in real time. The periodogram 
is simply the square of the magnitude of the Fast Fourier 
Transform (FFT) of a finite duration sequence x(n) . 

1. Discrete Fourier Transform 

Let x(n) be a discrete periodic signal of period N, 
then the definition for the discrete Fourier series pair is 
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JV-1 



xln] =52 aj^e 



jk{2%/N) n 



U) 



k‘0 



N-l 



3 x[n] 

k ^ 



( 2 ) 



n=Q 



It can be clearly seen in equation (1) that the discrete 
Fourier series is a finite length of N terms, compared to the 
infinite number of terms found in the continuous Fourier 
series. In addition, the largest number of terms present is 
N, hence {a^} will be a periodic sequence with period N. 

There is no discrete Fourier representation for 
nonperiodic discrete signals. Nonperiodic signals will 
normally be of infinite length and are represented by the 
following discrete-time Fourier transform pairs; 

Do*2it 

x[nl = — f Xia)e^“''da (3) 

L 

XiCl) = ^ x{n\ , ( 4 ) 

where fio is any real number. The discrete-time Fourier 
transform of x[n] is X(n) which is a continuous function 
having a period of 2n. 

Since x[n] is normally infinite in length and X(fl) is 
continuous, it is not practical to implement the discrete-time 
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Fourier transform on the computer. One approach to solving 
this problem is to make x[n] nonzero only for 0<n<W-l. This 
approach lead to the development of a new transform called the 
discrete Fourier transform (DFT) given by 

N-l 

x[n] = — 53;f[Jc] 0 ^ n s 2^-1 (5) 

N-l 

X[k] = 5^x[n] 0 i k ^ N-l , (6) 

n=o 

where X[kJ is the DFT of x[n], [Ref. 6:p.79]. 

2. Fast Fourier Transform 

There are a total of N terms involved in the summation 
of equation (6) and its computation requires (N-l) complex 
additions. In addition, each computation in the summation 
requires one complex multiplication. To compute a N-point DFT 
of a N-point sequence in a straight forward fashion, N(N-l) 
complex additions and multiplications are required. This 
is a significant computational strain for a computer and much 
effort was expended to develop ways of making the computation 
of the DFT faster. The early to mid 1960 's saw great progress 
in this area. The Cooley-Tukey algorithm, the Good-Thomas 
algorithm, Rader algorithm and chirp z algorithm (all which 
came to be known as Fast Fourier Transforms) , just to name a 
few, resulted in increased computational speed over the DFT. 
For a more complete listing of FFT algorithms as well as their 
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developments see [Ref. 6;pp. 90-115]. All of these algorithms 
offer features that are particularly useful and efficient when 
implemented on certain classes of computers. 

Regardless of how X(k) is computed, once it is 
obtained, the periodogram spectral estimate is computed via 

S^(Jc) = jX(Jc) = X*(k)X(k) , k = 0,1, . . .N-1, (7) 



where X'(k) is the complex conjugate of X(k) . 



C. BLACKMAN-TUKEY SPECTRAL ESTIMATION 

The Blackman-Tukey spectral estimator was named after R. 
Blackman and J. Tukey who published their work in 1958. The 
procedure estimates the autocorrelation function (ACF) as 






N-l-k 

— ^ k^O , 1 , . . . , iV-1 

i ^ /l-O 

f^E-k] k=-(i\T-l) ,-(^-2) , . . . ,-l 



( 8 ) 



which is Fourier transformed to obtain the PSD estimate. The 
periodogram is given by 

= E . (») 



The inferior performance of the ACF estimator is the cause for 
inferior performance of the Blackman-Tukey estimator. The ACF 
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estimator is biased regardless of sample size and has an 
expected value of 



E[f^[i:]] = \k\^N-l . (10) 



Examination of equation (10) demonstrates that the mean value 
of the ACF estimator is equal to the true ACF weighted by a 
Bartlett (triangular) window. An unbiased ACF could be used 
by replacing the 1/N factor in equation (8) with l/(N-l/c|). 
This could result in a negative spectral estimate since a 
positive semidefinite sequence can not be guaranteed when 
using an unbiased ACF estimator. This estimation approach was 
most popular until the development and implementation of the 
FFT algorithm. [Ref. 7]. 

D. WINDOWING EFFECTS 

Windowing (often referred to as weighting or shading) is 
an important topic that is integral to all classical 
estimation methods. We use windows to control the effects of 
sidelobes associated with these classical spectral estimators. 
"Tapering functions" is another name often given to data 
windows. The fundamental purpose of a data window is to 
lessen the bias in periodogram estimates. [Ref. 4;pp. 136- 
146] . 

Processing of finite duration data sets presents unique 
problems in analyzing the harmonic characteristics of that 
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data. It is important to pay particular attention to 
detecting spectral components in the close proximity of 
stronger spectral components as well as to resolvability of 
those components. The FFT assumes sequences to be periodic. 
That is the sampled data being analyzed is one complete period 
of an infinitely long periodic sequence. Of all possible 
frequencies in the data set only those that are located at an 
FFT bin center will be projected as a unique value in the 
frequency domain. All other frequencies will have non-zero 
projections over the entire frequency domain. This is termed 
spectral leakage and is a consequence of data records of 
finite duration [Refs. 4 and 5]. 

In nearly all cases of interest, the spectral components 
present in observed data will be at frequencies other than 
those located at FFT bin centers. Frequency components not 
located at a bin center will be non periodic in the 
observation window. This results in discontinuities at the 
observation boundaries, which in turn results in leakage over 
the entire range of the FFT. 

Consider a finite data record that is a portion of an 
infinite sequence that has been multiplied by some window. 
Consider an observed data sequence Xo[n] of N points that is 
the product of a rectangular window, rec(n) , of unit amplitude 
and an infinite-duration sequence, x(n) , where 
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f 1 73=0,1,2, . . . , (N-1) 

rec(n) = \ 

[ 0 otherwise 

and the observed data sequence is 

Xq [33] = x[n] • rec[n] 

The assumption made here is that all unobserved samples are 
essentially zero. Therefore, data that is processed "as is" 
is data that has been rectangularly windowed. Window 
functions other than the rectangular window function are 
called weighting functions that reduce the effects of boundary 
discontinuities. This smoothing of the observed data at the 
boundaries is one goal of windowing. 

Multiplication of the time series in equation (12) results 
in the convolution of the transfoinns of x(n) and rec(n) in the 
frequency domain. This results in a broadening or smearing of 
the power of the spectral components into adjacent frequency 
bins, provided the spectral component is not located at a bin 
center. If non-rectangular windows are used smearing will 
occur regardless of spectral component location. The 
narrowest spectral response of a windowed sequence can be no 
less than that deteirmined by the mainlobe of the transform of 
the window, independent of the data [Ref. 4:p.l37]. These 
mainlobe widths are different for each window function. For 
example, the main lobe width (half-power bandwidth, which is 
3 dB down from the peak response) for a rectangular window is 



( 11 ) 



( 12 ) 
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about the reciprocal of the observation interval (NT, where T 
is the time interval between samples) . The leakage mentioned 
earlier results in sidelobes of the main spectral components. 
These sidelobes may bias the amplitudes of adjacent frequency 
responses and may completely mask the presence of weaker 
signals, thus preventing their detection. 

Another goal of tapered windows is to achieve better 
sidelobe levels than those of the rectangular window. By 
decreasing the sidelobe levels the bias can be reduced. The 
drawback here is that sidelobe depression can only be 
accomplished by broadening the windows mainlobe frequency 
response. This in turn reduces the spectral resolution. 
Table I presents performance characteristics of different 
window functions. The trade-off between sidelobe depression 
and mainlobe resolution must be considered when determining 
the optimum window of choice. 

Table I. Window Performance Characteristics 

[Ref. 4: p. 143]. 



WINDOW 

NAMES 


HIGHEST 

SIDELOBE 

LEVEL 


1/ 2-POWER 1 
BW 1 

(DFT BINS) 1 


Rectangle 


-13.3 dB 


0.89 


Triangle 


-26.5 dB 


1.28 


Hann 


-31.5 dB 


1.44 


Hamming 


-43.0 dB 


1.30 
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E. ZERO PADDING 



There are two primary purposes for zero padding a data 
set. First, to allow for transform value interpolation 
between the original N transform values. To illustrate, let 
us consider a data set x(n) of a finite length N. Next, zero 
pad that data set with N zeroes such that 



[x(n) n=0,l,2, AN-1) 

yin) = < 

I 0 n=N, {N+D {2N-D . 



(13) 



From equation (6) the DFT of the 2N point data sequence y(n) 
is 

2N-1 

= Yr yin] e-J*< 2 «/ 2 Mn 

13*0 

N-1 

= Jc=0 , 1, 2 , . . . , (2i\T-l) . (14) 

Z3«0 



Comparing the results of equations (6) and (14) note that 



' X[k/2] ic=0,2,4. 



Ylk] 



i N-l 
J3*0 



ic=l,3,5. 



(2N-2) 

(2N-1) 



(15) 



and that the 2W-point DFT of y(n) is identical to the N~point 
DFT of x(n) at even index values and therefore the odd values 
of k represents the interpolated DFT values between the 
original W-point DFT [Ref 4:pp. 43-44]. The second purpose of 
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zero padding is to augment a data set to ensure it is of a 
length that is a power of two to allow the use of an FFT. 

To demonstrate the effects of zero padding consider a 16- 
point rectangular window. The DFT of a rectangular window 
will result in a digital sine function of the form 



D (f) 



' sinjnfTN) 
sininfT) 



(16) 



The digital sine function (D^) possesses sidelobes. 
Nevertheless, when we calculate and plot the transform we 
notice only a central spike at the origin with no sidelobes 
evident (Figure 1) . The reason that no sidelobes are visible 
is because the FFT of the non-zero-padded rectangular window 
series examines at its zero crossings and therefore, the 
structure of the sidelobes remains hidden. The zero crossings 
of occur at FFT bin centers and by zero padding one can 
calculate interpolated values between bin centers which allows 
the presentation of the sidelobe structure. Figure 2 
represents the magnitude of the FFT of a 16-point rectangular 
window zero-padded to 32-points. It is important to 
understand that zero-padding a data sequence prior to the DFT 
will not improve the resolution of the periodogram. This is 
the same principle that applies to more typical spectral 
estimation problems. 
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Figure 1. Magnitude of the FFT output of a 16-point 
Rectangular Window 
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bin centers 



FFT MAG OF 16-PT RECT WIN ZERO-PADDED TO 32-PT 




Figure 2. Magnitude of the FFT output of a 16-point 
Rectangular Window zero-padded to 32-point 
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frequency point 



To better illustrate the effects of zero-padding 
consider a sinusoid of unit amplitude intermixed in Gaussian 
white noise. The noise has a variance of 0.0005 which 
corresponds to SNR of 30 dB defined by 



SNR - lO-log^o 



2o2 



(17) 



where A = sinusoidal amplitude and is the variance of the 
Gaussian noise. In this illustration the data record length 
(N) is 32 and sampling frequency {/,) is 32 Hz. The 
sinusoidal frequency used is 10 Hz. The bin centers of th/^ 
FFT will fall at integer multiples of fj/N. In this example, 
fj/N = 32/32 = 1 Hz per bin. The sinusoidal frequency of 10 
Hz is located exactly at a FFT bin center. The spectral peak 
is clearly displayed (Figure 3) since the frequency is at bin 
center. Also note virtually no evidence of the presence of 
sidelobes since the digital sine function was interrogated at 
the zero crossings. Figure 4 is the results of the same 
signal described above with the exception of being zero padded 
to 64-points. Note that the sidelobes of Df, are now visible 
because of the zero-padding. 

Now let us examine the effect of a sinusoidal frequency 
not located at bin center. Let the sinusoidal frequency be 
10.7 Hz. Figure 5 is the result of a non-zero-padded 32-point 
sequence. The sidelobes are present even though the sequence 
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Figure 3. Periodogram of sinusoid in Gaussian noise, f=10 Hz 
(bin center) , 32-point sample, non-zero-padded, f=32 Hz 
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Figure 4 . Periodogram of sinusoid in Gaussian noise, f=10 Hz 
(bin center) , 32-pt sample zero-padded to 64-pt, f=22 Hz 
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Figure 5. PeriocJogram of sinusoid in Gaussian noise, 
Hz (not at bin center) , 32-pt sample non-zero-padded, 



f=10,7 
>3 2 Hz 
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was not zero-padded. This is because the sidelobes of Dj^ are 
now visible since the FFT is interrogating at locations 
other than its zero crossings. Spectral leakage is also 
evident, which has smeared signal energy into nearby frequency 
bins. Notice a broader, less-pronounced mainlobe results. 

F. STATISTICAL PROPERTIES OF THE PERIODOGRAM 

The distribution of a white Gaussian noise data set x(n) 
is given by: 

x{n) ~ . (18) 

Consider the case where the mean of x(n) is zero (i.e. ii=0) 
and the sample is of size W. The DFT of x (n) , defined as 
X(k) , will be comprised of real and imaginary parts, labeled 
A(k) and B(k) respectively, that are orthogonal linear 
combinations of x(n ) . For simplicity let X(k) be normalized 
by 1/SQRT(W) . Therefore, both A(k) and B(k) are Gaussian 

random variables with the distribution N(0,aj) and are 
mutually uncorrelated. The sum of the squares of the real and 
imaginary parts of X(k) is the periodogram (Per^) of x(n) and 
is given by 

Per^ik) =A^{k)+B^{k) . (19) 
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The probability density function of Per^(k) is chi- 
squared with 2 degrees of freedom. The mean and variance of 
Perj^(k) are 



E [Per^(k)] = 2o^ for all k 



Var \Perj^[k)] 



4ol; 


*-o,| 


00 

Q 


.=o.| 



( 20 ) 



( 21 ) 



For a derivation of equations (20) and (21) see 
[Ref. 8]. 



27 



III. SPECTRAL ESTIMATION USING THE KALMAN FILTER 
A. BACKGROUND 

As described in Chapter II the problem with the use of 
FFT-based spectral estimation is the trade-offs associated 
with sidelobe suppression versus spectral resolution. The 
application of tapered windows to a time series data set for 
the purpose of suppressing the sidelobes results in the 
broadening of the spectral mainlobe. Each tapered window 
possesses different performance characteristics as presented 
in Table I. It was demonstrated by Go [Ref. 5] that the 
application of the Kalman filter to the periodogram of an 
unwindowed spectral record for minimizing the variance with 
minimal degradation of the spectral resolution showed promise. 
Figure 6 allows the comparison of the resolution maintained by 
the Kalman filter relative to the mainlobe broadening seen in 
the Hamming window. The conclusions reached by Go were: 

1. The Kalman filter can enhance the spectral peaks of the 
periodogram of an unwindowed time series. 

2 . The resolution for the multi-spectral peaks of the 
periodogram is largely preserved by the Kalman filter. 

3. Reliable signal detection was achieved for signals with 
SNR's (after processing) of 12 dB for ^ in the range 
between 100,000 and 700,000. 

4. Signal detection was possible for SNR's (after 
processing) down to 6 dB in a small percentage of noise 
realizations. 



28 







o 

K) 



O 

CM 



d* 

c 

0 

3 

cr 

t! 







o 

lO 



o 

cs 



Qp opn)(u5Duj 



gp opn)|u5ouj 



Figure 6. Periodogram of f=10.7 Hz, input SNR=20 dB, 
filter using Hamming window and Kalman filter jS=80K 
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Kalman filtering is an estimation or prediction procedure 
which was introduced by R.E. Kalman and R.S. Bucy in the early 
1960 's. The basic idea of a Kalman filter is to recursively 
update the estimate of the state of a system by comparing 
these estimates to the system measurements. A new estimate 
follows each measurement. The basic ingredients of the Kalman 
filter are a stochastic differential system (called the system 
model) and a set of observations (called the measurement 
model) which allow for the computation of the state vector 
estimates. The measurement and the state estimate need not be 
of the same dimensionality. [Ref. 9]. 

The Kalman filter algorithm demonstrated by Go was written 
by Dr. Roberto Cristi at the Naval Postgraduate School, 
Monterey California in 1988. The algorithm was developed to 
detect piecewise constant segments of time series data 
corrupted by noise. If data under examination consists of 
piecewise constant segments, then there are two possibilities 
at each new observation. They are; 

1. the current observation is a continuation of the previous 
piecewise constant segment of data or 

2. the current observation is the first element of a new 
segment of data. 

The algorithm uses a filter parameter )3 to quantify the 
likelihood of a transition from one piecewise constant segment 
to another. If no transition is detected then the current 
observation is filtered using a Kalman filter updated with the 
current Kalman gain. If a transition is detected then the 
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current observation is filtered using a reinitialized Kalman 
filter. The sensitivity of the filter to detect transitions 
between piecewise constant segments is a function of the 
parameter If the value of /3 is "too small", then the 

filter will reinitialize too often resulting in less smoothing 
of the data. If the value of ^ is "too large", then over- 
smoothing will occur and transition points will not be 
detected. For a detailed description of the development of 
this algorithm see [Ref. 5: ch. 3]. To demonstrate the 

effects of different values of j8 consider a sinusoid corrupted 
by Gaussian noise. Figure 7 shows the result of filtering a 
periodogram of an unwindowed sinusoid (freq=10.7 Hz) time 
series with a sample size of 128-pts zero padded to 256-pts 
sampled with (sampling frequency) =64 Hz and an input SNR=-9 
dB. Values for ]8 of 8k, 80k and 800k (where k=10^) where used 
to filter this data. Note the added detail in the plot of the 
filter output for /3=8k versus the over smoothing where jS=800k. 

B. EXPERIMENTAL DATA 

The data used in this research was computer generated. 
It is sinusoidal data in additive Gaussian white noise of 
different variances. Every data set had at least one sinusoid 
of unit amplitude with a preselected SNR. The noise variance 
(a^) for that data set was calculated using equation (17) . 
For data sets with more than one frequency component, the 
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Figure 7. Periodogram of 10.7 Hz sinusoid of -9 dB input SNR 
filtered for jS values of 8k, 80k, 800k respectively. 
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amplitude of each component was the determining factor in 
achieving that component's desired SNR. 

Transform lengths of 128-pts zero-padded to 256-pts were 
used exclusively in this research. Zero-padding does not 
improve FFT processing gain, which for best case can be 
approximated by 

[log 2 (data record length) - 1] x 3 dB (22) 



provided the data is stationary over the data length. The 
processing gain for an 128-pt FFT is about 18 dB, therefore 
the realized SNR on the output side of the FFT is 



+ PG 

= + 18 



(23) 



where SNRg^ is the SNR on the output side of the FFT and SNRj„ 
is the SNR of the input signal. Data used in this study will 
be of this form unless otherwise specified. 

C. KALMAN FILTER MODIFICATIONS 
1. Dynzunic Beta 

It was shown [Ref. 5] that the Kalman filter algorithm 
can be successfully adapted for smoothing spectral data. The 
filter parameter /? is the controlling factor that determines 
how much the noise portions of the periodogram are smoothed 
and how well the transition points are detected. The result 
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for an appropriately chosen /3 will be a smoothed periodogram 
with narrow mainlobes that are essentially unaltered from the 
original periodogram. A common trait of the original Kalman 
filter is a tapering (i.e., roll-off) effect on the higher 
frequency side of the mainlobe. This is a result of the 
Kalman filter failing to reinitialize after the peak and on 
the down-transition of the periodogram. A smaller valued /3 
would be more sensitive and would more closely follow the 
down-transition of the original periodogram. Unfortunately, 
as shown in Figure 7, small values of /3 provide less smoothing 
and little improvement over the original periodogram. 

This is the motivation for allowing to be dynamic 
during the filtering process. An optimum value of /3 is used 
during the initial filtering process until a transition is 
detected. Once a transition is detected the value of )3 is 
decreased to insure maximum sensitivity ()8=1 is used in this 
study) , which minimizes the tapering effect on the down- 
transition side of the spectral component. 

Consider a sinusoid where f=10.7 Hz and SNRj„=-9 dB. 
The periodogram of the unwindowed data record is filtered 
using the original Kalman filter algorithm with /3=80k. The 
same periodogram is again filtered using the Kalman filter 
with the dynamic /3 (i.e., /3=80k or 1). Figure 8 shows the 
result of that type of processing. Notice the tapered 
response on the higher frequency side of the spectral peak in 
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Figure 8. Periodogram of 10.7 Hz sinusoid of SNRi„=- 
filtered with constant /3=80k and a dynamic /3=80k 
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the case for constant )8 filtering. The shape of the 
unwindowed periodogram on the high frequency side of the 
mainlobe is clearly more preserved in the dynamic /3 case. 
Another important observation is that there is no change in 
the dB level of the noise floor for either filter output. 

It is not uncommon for the dynamic /3 filter to provide 
similar results as when using the constant /3 filter for some 
noise realizations. While /3 is at the "very" sensitive 
setting and the data does not exceed a given transition 
threshold then the algorithm will redefine 0 to the original 
filtering value. This redefining of may occur in the 
vicinity of the spectral peak. This can result in filtering 
of the down side of a spectral peak with the same value of /3 
used in the constant jS filter. 

Although there are instances when the dynamic /3 will 
not improve filter performance, it was exceedingly rare (less 
than 10 out of 256 trials) for the dynamic filter to provide 
less satisfactory results than the constant /3 filter. The 
LOFAR output is the presentation of 128 single step traces. 
Even if only a few of the 128 traces are improved by the 
dynamic filter process, then those improvements could lead to 
enhancement of signal detectability and resolution. 

2 . Data Svibstitution 

The second modification made to the Kalman filter 
algorithm was driven by similar motivation as for the dynamic 
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jS filter discussed earlier. It is clear that the Kalman 
filter can smooth spectral data while leaving the mainlobe of 
the periodogram virtually unaltered for appropriately chosen 
/3 and SNR^, of 9 dB and larger. The second modification was 
made in the attempt to improve upon the filters ability to 
preserve resolution of the original spectral data. In the 
frequency domain the spectral component of a sinusoid will 
span only a small number of FFT bins as compared to the total 
number of bins/ frequencies represented in the FFT output. For 
an appropriately chosen /? the total number of times the Kalman 
filter detects transitions from one piecewise smooth segment 
to another is few compared to the number of times the filter 
detects no transitions. 



series and let X(k) be the output of the Kalman filter (for 
either the constant j8 or dynamic j8 algorithm) computed from 
the filtering of Y(k). The Kalman filter was modified to 
produce a new output (Xnf(k) ) 



This modification can be implemented in conjunction with 
either the dynamic or constant /3 filters. The notation 
"modified” Kalman will be used to describe this modification 
throughout this study. 



Let Y(k) be the periodogram of an unwindowed time 




no transition detected 



up- transition detected . (24) 
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Figure 9 shows the Kalman output of a dynamic /3 
process versus the output from the dynamic /3 algorithm 
modified as described above in equation (24) . Notice there is 
not an obvious improvement in resolution of the spectral peak 
over the dynamic only case (upper right Figure 9 ) , but there 
is a measurable lowering of the noise floor realized for the 
dynamic /3 filter modified to substitute raw data in for the 
estimate during detected transitions. The noise realization 
filtered in Figure 9 displays about a 1.5 dB improvement in 
raising the spectral peak height above the noise floor. This 
is a consistent result observed throughout this study. 

D. KALMAN VERSUS HAMMING 

Results presented by Go [Ref. 4] indicated that the Kalman 
filter appeared to outperform the Hamming window in signal 
detection at SNR^ of 9 dB and below for /3=300k. The 
relatively small number of noise realizations looked at in 
reference 4 were insufficient to estimate just how much better 
Kalman was over the Hamming window in signal detection. Here 
the outputs of the Hamming window and Kalman filters (for 
/3=80k,135k and 300k) of both the dynamic /3 and the "modified” 
Kalman filter using dynamic jS were examined for 128 noise 
realizations. This was done twice. Once for SNJ?i,=-9 dB and 
once for SNRj„=~12 dB. For both cases the sinusoidal frequency 
was f =2 0.0625 Hz. The output for each noise realization was 
examined to determine a hit or miss for signal detection for 
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Figure 9. A comparison of the performance of Kalman filtering 
(dynamic /S vs modified dynamic /3) of an periodogram of an 
unwindowed time series of f=10.7 Hz, SWi?^=-9 dB 
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each filter process. An example of one noise realization for 
SNR^=-9 dB for the Hamming window (upper right) , Kalman with 
dynamic ^=80k (lower left) , and Kalman "modified” filter with 
dynamic /3=80k (lower right) is provide with Figure 10. If the 
output displayed a peak at the true spectral location then a 
detection was recorded for that filter process and that noise 
realization. The totals were tabulated and are presented in 
Table II. 

Table II. Detection comparisons for Kalman vs Hamming for 
128 noise realizations at SWi?^=-9 dB and for 128 noise 
realiza-tions at SNRj„=-12 dB. The number of misses are 
provided . 



1 Filter Process 


SWJ?i,=-9 dB 


SNR^=-12 dB 1 


Hamming Window 


31 


59 


Dyn /3 = 80k 


19 


52 


Mod Kal, Dyn jS=80k 


19 


49 


Dyn /S = 135k 


24 


55 


Mod Kal, Dyn /5=135k 


24 


52 


Dyn |8 = 300k 


36 


70 


Mod Kal, Dyn /S=300k 


36 


67 



A statistic that was not collected was the number of false 
detection (i.e., false alarms) made by each different 
filtering process. Since the spectral location of the false 
alarms are random and independent of the noise realization, 
and each trace is one of 128 that will be displayed in a LOFAR 
output, the false detections should have negligible impact. 
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Figure 10. Periodogram (f=20.0625 Hz, SNRj„=-9 dB) filtered 
using Hamming window, Kalman Filter Dyn /3=80k and Modified 
Kalman Filter Dyn jS=80k 
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E. LOFAR OUTPUT 

As mentioned earlier, for very long signal segments it can 
be cumbersome and time consuming to look at the periodograms 
or filtered periodograms trace by trace. Also at SNR^ of 6 
dB or less it is sometimes difficult to distinguish signal 
from noise for many of the individual traces. Another problem 
of examining spectral data trace by trace is the difficulty of 
detecting the dynamics that may be present in some signals. 

A solution to these problems is to display all (or at 
least a substantial number) of the periodograms simul- 
taneously. This display needs to be relatively fast and 
compatible with the SUN workstations. The procedure used is 
to stack the periodograms into a matrix that is 128 by 128 in 
size. The largest numerical value in the matrix is assigned 
the darkest grey-tone and conversely, the lightest grey-tone 
is assigned to the lowest value. Sixteen intensity levels are 
used for the mapping of the remaining data. To each intensity 
level is a corresponding equal range in the spectral density. 

The graphics package used to create this image is NCAR. 
NCAR Graphics is a collection of FORTRAN 77 programs and 
subroutines that can be used to generate and plot computer 
graphics suitable for the display of scientific data. NCAR 
Graphics conforms to the Graphical Kernel System (GKS) 
standard. Level OA (zero A). 
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The LOFAR outputs used in this study are created using the 
subroutine HALFTON which draws a halftone picture from data 
stored in square or rectangular arrays. A time domain 
(sinusoidal) signal corrupted by additive Gaussian white noise 
is generated. The signal length is 16384 points. The signal 
is processed in segments using 128 non-overlapping contiguous 
samples, zero-padding to 256 points. Each segment provides 
one periodogram which is stacked in array form. The resulting 
array is a square matrix 128 by 128, where the rows are the 
periodograms of the individual segments and the columns 
represent frequency. The matrix is read into HALFTON and is 
displayed as a LOFAR output. Figure 11 is an example of a 
LOFAR output of rectangular windowed data with f=10.7 Hz, 
SWR^=-12 dB and /,=64 Hz. The horizontal axis is the frequency 
axis which goes from 0 to 32 (/,/2) Hz while the vertical axis 
represents time running from the bottom to the top. The 10.7 
Hz line is quite obvious even though the SNR^ is only 6 dB. 
This dramatic increase in detectability over the individual 
periodograms is due to the additional processing gain obtained 
by the visual display. 
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LOFAR gram. RECT window. f=10.7 Hz. 5NRout=6 dB 




Figure 11. LOFAR output of rectangular windowed signal 
f=10.7 Hz, SNR^=-12 dB and f,=64 Hz 
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IV. EXPERIMENTAL RESULTS 



A. EFFECTS OF FILTER MODIFICATIONS 
1. Dynamic versus Constant 

Allowing /3 to be dynamic during the filtering process 
often results in a more accurate representation of the 
spectral shape on the higher frequency side of the mainlobe 
than achieved by the constant filter. Figure 12 thru 
Figure 21 are 10 different noise realizations which compare 
the performance of a constant /3 filter with a dynamic /3 
filter. The spectral data has two frequency components. One 
is fl at 10.7 Hz at an SNRj„ of -9 dB and the other one is f2 
equal to 11.7 Hz at an SNRi„ of -12 dB. A 1 Hz separation is 
equivalent to 4 bin separation for a sample size of W=128 at 
a sampling frequency of 64 Hz. The dynamic /S filter provided 
more accurate representation of the original periodogram in 4 
of the 10 noise realizations without a single case of inferior 
performance. 

Figure 22 shows two LOFAR outputs allowing the 
comparison of the outputs for constant and dynamic /3. In this 
representation it is not clear that there is any 
distinguishable difference between the two filtering outputs. 
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Figure 12. Periodogram, Noise Realization 1 (f 1=10. 7 Hz at 
SWR^=-9 dB and f 2=11. 7 Hz at SWR^=-12 dB) . Kalman filter using 
a constant )8 and a dynamic j8. 
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Figure 13. Periodogram, Noise Realization 2 (f 1=10. 7 Hz at 
SNRi„=-9 dB and f2=11.7 Hz at SNRj„=-12 dB) . Kalman filter using 
a constant /3 and a dynamic /S. 
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Figure 14. Periodogram, Noise Realization 3 (f 1=10. 7 Hz at 
SNRi„=-9 dB and f 2=11. 7 Hz at SNRi„=-12 dB) . Kalman filter using 
a constant j8 and a dynamic /3. 
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Figure 15. Periodogram, Noise Realization 4 (f 1=10. 7 Hz at 
SNRi„=-S dB and £2=11.7 Hz at SWRj„=-12 dB) . Kalman filter using 
a constant /3 and a dynamic /S. 
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Figure 16. Periodogram, Noise Realization 5 (f 1=10. 7 Hz at 
SNR^=-9 dB and £2=11.7 Hz at SNRj„=’~12 dB) . Kalman filter using 
a constant and a dynamic 
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Figure 17. Periodogram, Noise Realization 6 (f 1=10. 7 Hz at 
SNRj„=~9 dB and f2=11.7 Hz at SNRi„=~12 dB) . Kalman filter using 
a constant j8 and a dynamic /3. 
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Figure 18. Periodogram, Noise Realization 7 (f 1=10. 7 Hz at 
SNRi=-9 dB and £2=11.7 Hz at SNRi„=^-12 dB) . Kalman filter using 
a constant /3 and a dynamic /3. 
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Figure 19. Periodogram, Noise Realization 8 (f 1=10. 7 Hz at 
SNRj„=-9 dB and f2=11.7 Hz at SNRj„=-12 dB) . Kalman filter using 
a constant jS and a dynamic /3. 
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Figure 20. Periodogram, Noise Realization 9 (f 1=10. 7 Hz at 
SNRj„=-9 dB and f 2=11. 7 Hz at SNRj„--12 dB) . Kalman filter using 
a constant and a dynamic /3. 
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Figure 21 . Periodogram, Noise Realization 10 (f 1=10. 7 Hz at 
SNRj„=-9 dB and f2=11.7 Hz at SNi?^=-12 dB) . Kalman filter using 
a constant /3 and a dynamic j8. 
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Figure 22. LOFAR presentation of a two frequency component 
signal (f 1=10. 7 Hz at SNRi„=-9 dB and f 2=11. 7 Hz at SNRi=-l2 
dB) . Kalman filter output using a constant and a dynamic /3. 
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2. Kalman versus "Modified" Kalman 



Modifying the Kalman filter to substitute raw spectral 
data in place of filter estimates after up transitions are 
detected provides two important benefits. First, in single 
trace outputs more detail in the vicinity of the spectral 
mainlobe is preserved in a large number of the noise 
realizations. Second, the noise floor is consistently lowered 
anywhere from 1 to 3 dB, with 1.5 dB appearing to be about the 
average improvement. This result can have a dramatic effect 
on the quality of the LOFAR image output. Figure 23 thru 
Figure 32 are 10 different noise realizations of periodograms 
filtered using the dynamic /3 and the "modified" dynamic ^ 
algorithms. The input time series contains a sinusoid at a 
frequency f at 10.7 Hz with a SNRj„ of -9 dB. In 6 of the 10 
noise realizations is a measurable improvement when using the 
"modified" dynamic /S filter in that the distance of the 
spectral peak separation from the noise floor increases. In 
addition there are 3 noise realizations where the spectral 
resolution is better preserved by the "modified" dynamic /8 
algorithm than by non-modified dynamic 0 filter algorithm. 

Figure 33 shows the LOFAR output for the "modified" 
filter and the non-modified filter. Both filters used a 
dynamic /S of 135000. The reduction of the noise background is 
obvious. The "modified" dynamic filtering scheme shows 
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promise to improve signal detectability as well as the 
frequency resolution. 
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Figure 23. Periodogram, Noise Realization 1, (f 1=10. 7 Hz at 
SNRi„=-9 dB) . Kalman filter performance using a dynamic /8 and 
"modified" dynamic /S of 135000. 
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Figure 24 . Periodogram, Noise Realization 2, (f 1=10. 7 Hz at 
SNRj„='-9 dB) . Kalman filter performance using a dynamic /3 and 
'•modified” dynamic /S of 135000. 
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Figure 25. Periodogram, Noise Realization 3, (f 1=10. 7 Hz at 
5WR^=-9 dB) . Kalman filter performance using a dynamic and 
••modified” dynamic jS of 135000. 
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Figure 26. Periodogram, Noise Realization 4, (f 1=10. 7 Hz at 
SNRj„=-9 dB) . Kalman filter performance using a dynamic jS and 
"modified" dynamic /8 of 135000. 
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Figure 27. Periodogram, Noise Realization 5, (fl=10.7 Hz at 
SWR^=-9 dB) . Kalman filter performance using a dynamic /3 and 
"modified” dynamic /S of 135000. 
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Figure 28. Periodogram, Noise Realization 6, (f 1=10. 7 Hz at 
5NJ?j,=-9 dB) . Kalman filter performance using a dynamic jS and 
"modified" dynamic |5 of 135000. 
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Figure 29. Periodogram, Noise Realization 7, (f 1=10. 7 Hz at 
SNRj„=-9 dB) . Kalman filter performance using a dynamic /3 and 
"modified" dynamic jS of 135000. 
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Figure 30. Periodogram, Noise Realization 8, (f 1=10. 7 Hz at 
SNRj„=-9 dB) . Kalman filter performance using a dynamic j8 and 
"modified” dynamic jS of 135000. 
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Figure 31. Periodogram, Noise Realization 9, (f 1=10. 7 Hz at 
SNRj„=-9 dB) . Kalman filter performance using a dynamic ^ and 
•'modified” dynamic /3 of 135000. 
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Figure 32. Periodogram, Noise Realization 10, (f 1=10. 7 Hz at 
SNRi„=-9 dB) . Kalman filter performance using a dynamic /8 and 
"modified” dynamic |8 of 135000. 
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Figure 33. LOFAR output (fl=10.7 Hz at SNRj„=~9 dB) . Kalman 
filter output using a dynamic and a "modified" dynamic /3 of 
135000. 



69 



B. "MODIFIED" KALMAM FILTERING VERSUS HAMMING WINDOWING 

It has been demonstrated that the Kalman filter more 
accurately preserves the spectral resolution of the original 
periodogram whereas the Hamming window will broaden the 
mainlobe (see Figure 6) . The statistics tabulated in Table II 
(p. 40) clearly show that the Kalman filter (/3=80k and 135k 
for both the “modified” dynamic /3 and the dynamic /3) is more 
successful at signal detection at an SNRj„ of -9 and of -12 dB 
than the Hamming window. 

Figure 34 displays the LOFAR output of a 3 component 
sinusoidal signal (f 1=10. 7 Hz at SWJ?^=-9 dB, f 2=14. 45 Hz at 
SNRj„=-15 dB and f 3=20. 25 Hz at SNRj„=~9 dB) of rectangularly 
windowed segments. Figure 35 permits the comparison of using 
the Hamming window and the "modified" Kalman filter with a 
dynamic of 135k on the same data. 

The "modified" Kalman filter preserved the spectral 
resolution of the original unwindowed spectra where as the 
mainlobe broadening is evident when applying the Hamming 
window. This broadening can have a major impact in detecting 
nearby spectral components. See Appendix B for a closer look 
at resolution comparisons between Hamming window and the 
"modified" Kalman filter. Also note the reduction in the 
noise background of the Kalman LOFAR output which enhances the 
-15 dB component more clearly than the Hamming window. 
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LOFAR GRAM. RECTANGULAR WINDOW 




Figure 34. LOFAR output of Rectangular windowed data of 3 
component signal (fl=10.7 Hz at 5WR^=-9 dB, f2=l4.45 Hz at 
5NRi,=-15 dB and f3=20.25 Hz at SNR^=-12 dB) . 
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Figure 35. LOFAR output of the 3 component signal of 
Figure 34. Hamming window and Kalman ("modified" dynamic 
/5=135k) filter. 
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C. EFFECTS OF SNR^ 

Kalman filtering consistently permits detection of signals 
corrupted by additive Gausian white noise at a SNR^^ of 6 dB 
or greater. Consider the two extreme cases of signal only and 
noise only. A signal at SNRj„=20 dB (signal only case) 
filtered using the "modified” Kalman filter (using a dynamic 
/8=8k, 80k and 800k) produces excellent results as shown in 
Figure 36 where the "modified" Kalman filter detects the 
frequency component and preserves the resolution regardless of 
the used. In the case of noise only (SWi?^=-60 dB) the 
"modified" Kalman filter output provided in Figure 37 is 
identical for all three jS's used in processing the original 
periodograms. 

To take a closer look at the "modified" Kalman filter 
performance with SNR^ of less than 9 dB, we consider a signal 
of 6 non-equally spaced frequency components at fl=3.3 Hz with 
SNRi=-9 dB, f2=8.7 Hz with SNRi„=-12 dB, f 3=14. 4 Hz with SNJ?^=- 
13 dB, f 4=19. 3 Hz with SNR^^-IA dB, f 5=23. 7 Hz with SNRi„=-15 
dB and f 6=28. 3 Hz with SNRj„=-16 dB. We filter this signal 
with a Haitiming window and a "modified" Kalman filter using a 
dynamic jS=80k. Figure 38 displays the LOFAR outputs of these 
filters. In both LOFAR outputs fl, f2 and f3 are detectable. 
However, note the reduced noise background in the "modified" 
Kalman filter and the superior resolution of the "modified" 
Kalman filter (most obvious for fl) . The tonal line of f4 
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Figure 36. Periodogram of f=10.7 Hz at SNRj„-30 dB 
only) filtered using a Kalman (''modified” dynamic /5) 
for j8=8k, /3=80k and |8=800k. 
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Figure 37. Periodogram of f=10.7 Hz at SWH^=-60 dB 
only) filtered using a Kalman ("modified” dynamic /8) 
for jS=8k, ]8=80k and /8=800k. 
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Figure 38. LOFAR outputs using Hamming window and "modified" 
Kalman (dynamic /3=135k) . Freqs(Hz) at SNRSj„(<iB) ; fl=3.3 at -9, 
f2=8.7 at -12, f 3=14. 4 at -13, f 4=19. 3 at -14, f 5=23. 7 at -15 
and £6=28.3 at -16. 
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{SNRj„=-14 dB) is also present in both but is not nearly as 
distinct. It is much more difficult to observe f5 or f6 using 
either filtering process. 

Further processing gain (for stationary frequency 
components) can be realized by averaging the columns of the 
periodogram arrays. The procedure used in this study is to 
sum the columns and divide the sum by the number of rows (128 
in our case). Let Y(j,k) be an input array of size J by W 
then X(k) is calculated as follows; 



SumY(k) 

SumYnormik) 

X{k) 



j 



-T. 

7«1 



Y{j.k) 

J 



SumYik) 



max (SumY(k) ) 

= 10 ‘log {SumYnorm (k) ) 



k=l , 2 , .... ,N. (25) 



Figure 39 is the output of this process for the rectangular 
window, Hamming window, Kalman (•'modified” dynamic /3=80k) and 
Kalman ("modified” dynamic /3=135k) . All six frequency 
components are clearly present in all the displays. It 
appears however that Kalman is superior in smoothing the noise 
variance between components. 

To take this experiment one step further a new signal was 
generated using the same frequencies as before but at lower 
SNR's (fl at SNRi=-16 dB, f2 at SNR^=-17 dB, f3 at SNR^=-18 dB, 
f4 at SNRf„=-19 dB, f5 at SNR^=-20 dB and f6 at SNR^=-21 dB) . 
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Figure 39. Averaging of LOFAR outputs using Hamming window 
and "modified" Kalman (dynamic /3=80k and 135k). Freqs(Hz) at 
SNRs i„ (dB) ; 3.3 at -9, 8.7 at -12, 14.4 at -13, 19.3 at -14, 
23.7 at -15 and 28.3 at -16. 
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For completeness Figure 40 shows the LOFAR outputs for the 
Hamming window and the Kalman ("modified" dynamic )3=13 5k) 
filters for this experiment. Neither process appears to have 
much success at detecting any of the six frequency components. 
Figure 41 is the output of the process defined in equation 
(25) . The "modified" Kalman filter is much less successful at 
smoothing the noise variance between frequency components in 
this SNR range. This is due to the tendency of the Kalman 
filter to track the noise at small SNR's (also see Figure 37). 

The tendency of the Kalman filter to track the noise at 
small SNR's is a result of the statistical properties of the 
FFT of additive white Gausian noise (WGN) . If the noise input 
(time series) is a white Gaussian sequence with zero mean and 
variance , then the magnitude squared of the FFT \Y (k) p will 
have a chi-squared distribution with 2 degrees of freedom with 
a variance o\ where the mean equals the standard deviation. 
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Figure 40. LOFAR outputs of Hamming window and "modified” 
Kalman (dyn /S=135k) filter. Freqs at fl=3.3Hz, f2=8.7Hz, 
f3=14.4Hz, f4=19.3Hz, f5=23.7Hz and f6=28.3Hz with SNRi„'s of - 
16, -17, -18, -19, -20 and -21 dB respectively. 
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Figure 41. Averaged LOFAR outputs. Hamming window, "modified" 
Kalman using dynamic jS of 80k and 135k. Freqs; fl=3.3Hz, 
f2=8.7Hz, f3=14.4Hz, f4=19.3Hz, f5=23.7Hz and f6=28.3Hz with 
of -16, -17, -18, -19, -20 and -21 dB. 
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V. SOMMARY 



A. CONCLUSIONS 

The Kalman filters examined in this study (constant jS, 
dynamic /3 and ''modified” filter) all demonstrated the ability 
to enhance the spectral peaks of a periodogram of an 
unwindowed time series. The dynamic /3 and "modified” dynamic 
j8 filters provide the best resolution preservation of the 
mainlobe with negligible broadening. The Kalman filter 
significantly outperformed the Hamming window at SNR's^ of 6 
and 9 dB for j8=135k and jS=80k. The false detections for these 
filters cause negligible degradation to the output when 
displayed as LOFAR gram. 

The selection of the filter parameter /3 is dependent on 
SNR^. In addition, the selection of the filter parameter /3 
can be adjusted to maximize filter performance if prior 
knowledge of the desired filter output is available. Using 
the "modified” filter with /3=135k results in reliable signal 
detection at SNR^ of 9 dB in the periodogram with minimal 
false detections. For |8=80k, reliable detections are also 
made at SNRg^ of 9 dB but at the cost of an increase in false 
detections. Prior knowledge of the composition of the 
expected signal of interest will have a major impact on the 
selection of |8. The larger the expected SNR^, the larger 0 
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can be to provide good smoothing with maximal signal 
detectability. For SNR^ of 6 dB, a /3 of 13 5k provide the best 
results. For SNR^ less than 6 dB, a /3 of 80k is the optimum 
choice while for SNR^ greater than 9 dB, a /3 of 3 00k is 
appropriate. 

Viewing Kalman filter outputs on LOFAR outputs has 
resulted in detection of signals with SNR^ down to 3 dB. The 
LOFAR output also demonstrates the capability for the Kalman 
filter to process time varying spectra. See Appendix C for a 
more detailed description. 

The "modified" Kalman filter (for either constant or 
dynamic jS) can enhance the spectral peaks of a periodogram by 
consistently reducing the noise floor by 1 to 3 dB over the 
non-modified filter. This reduction in noise can have a 
dramatic impact when displayed as a LOFAR output. 

B. FUTURE WORK 

The Kalman filter is effective at preserving the spectral 
resolution of the original periodogram. Modifying the filter 
algorithm to substitute the original periodogram data for the 
filter estimates during up transitions enhances the Kalman 
filter performance in signal detection, noise reduction and 
spectral preservation. With the recent interest in signals 
with narrowband mainlobe structure (called "SWATH" signals) , 
further modifications of the Kalman algorithm for specific 
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application to the processing of these "SWATH” components is 
an area of future work. 
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APPENDIX A. USERS GUIDE / COMPUTER CODE 



A. SIMULATION USERS GUIDE 

The Kalman filtering algorithms used in this study are 
written in Pro-Matlab (version 3.5h). All algorithms used to 
support the simulations are also written in Matlab with the 
exception of the grey-tone shading program which is written in 
FORTRAN 77. The simulations were run on a NPS Sun work 
station. The NPS Sun stations use the UNIX operating system 
via an Internet Networking System. 

A complete simulation (from data generation to LOFAR 
output) can be executed providing all of the following code is 
copied into the directory in which Matlab is started in. This 
simulation is limited to a single stationary frequency 
component signal processed using the "modified" Kalman filter 
with a dynamic beta. The following steps should be followed 
to run this simulation: 

1. Change to the directory that all applicable programs are 
residing in. 

2. Start matlab by typing 'matlabl'. 

3. At the matlab prompt type 'kfilter' . 

4. Follow the instructions on the screen. 

The entire simulation will take approximately 15 minutes for 
a 128 by 128 Kalman filtered array output. If a hard copy 
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output is desired follow the procedures provided in Appendix 



D. The program KFILTER.M can be modified to generate more 
complicated data for processing. 



B . COMPUTER CODE 

All the computer code needed to reproduce the work in this 
thesis is provided below. 

1. Kalman Filter Algorithms 

a. "Modified" Kalman with Dynamic fi 
KALMODl.M 

This program filters the array of desired periodograms 
of time series data using a "modified" kalman filter with 
a dynamic Beta. The entering arguments are; 

- <file name> of PSD array 

- sample size 

- sigma (already calculated if data was generated 
using "KFILTER.M" program) . 

- total number of (n»m) points in the array (already 
calculated if "KFILTER.M" was used) . 

- the kalman filter parameter "Beta" 

The subroutine outputs are; 

M= periodogram of filter output in dB 
C= periodogram of filter output in power 

% Define the function "KALMODl" 

function [ M , C ] =ka Imod 1 ( psd , samps ize , s igma , totpts , betah ) 



% KALMAN FILTER 

x=zeros (2 , sampsize) ; 
pointer=zeros (2 , sampsize) ; 

beta=betah; 
betal = 1 

sv2=sigma^2 ; 

% create loop for the number of rows in psd 
r ows=t otpt s / samps ize; 
for k=l:rows 

y=psd(k, :) ; 
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el=0. 0; 
e2=0. 0; 

dl=0. 0; 
d2=0. 0; 

x(l,l)=y(l) ; 
x(2,l)=y(l) ; 
tau=l . 0; 

% MAIN LOOP 

for t=l: (sampsize-1) ; 

kl=1.0/ (tau+1.0) ; 
xll=x(l,t)+kl*(y(t)-x(l,t) ) ; 
dll=dl-beta; 

ell=el+(1.0/(2.0*sv2) ) * ( (y (t+1) -xll) ^ 2 ) ; 
cll=ell+dll; 

k2=1.0; 

xl2=x(l,t)+k2*(y(t) -x(l,t) ) ; 
dl2=dl+beta; 

el2=el+ (1.0/ (2 . 0*sv2) ) * ( (y (t+1) -xl2) ^2) ; 
Cl2=el2+dl2; 

kl=0 . 5 ; 

x21=x(2,t) +kl*(y (t) -x(2,t) ) ; 
d21=d2-beta; 

e21 = e2+(1.0/(2.0*sv2))*((y(t+l)-x21)"2); 

C21=e21+d21 ; 

k2=l. 0; 

x22=x(2,t)+k2*(y(t)-x(2,t) ) ; 
d22=d2+beta ; 

e22=e2+(1.0/(2.0*sv2) )*( (y(t+l)-x22) ^ 2 ) ; 
C22=e22+d22 ; 

% UPDATE STATES IN DYNAMIC PROGRAM. 

if clKc21 

x(l,t+l)=xll; 
el=ell ; 
dl=dll ; 
cl=el+dl ; 

pointer (1, t+1) =1 ; 

tau = tau+l; 
beta=betah ; 

else 

x(l,t+l)=y(t) ; 
el=e21 ; 
dl=d21; 
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% 

% 



cl=el+dl; 
tau=2 ; 

pointer ( 1 , t+1) =2 ; 
beta=betal ; 

end 



end 



if c22<cl2 

x(2,t+l)=y (t) ; 
e2=e22 ; 
d2=d22 ; 
c2=e2+d2 ; 
pointer ( 2 , t+1 ) =2 ; 
beta=betal ; 

else 

x(2,t+l)=xl2; 
e2=el2 ; 
d2=dl2 ; 
c2=e2+d2 ; 
pointer ( 2 , t+1 ) =1 ; 
beta=betah ; 

end 



END MAIN LOOP 



BACKWARDS SMOOTHING AND SUBSTITUTION 



tau=l . 0 ; 



if cKc2 

out=l ; 

else 



out=2 ; 

end 

y (sampsize)=x(out, sampsize) ; 

for t=sampsize: -1: 2 

out=pointer (out, t) ; 
xout=x(out, t-1) ; 

if out==2 

tau=l . 0 ; 
y (t-l)=xout; 

else 

tau=tau+l; 
y (t-l)=xout; 

end 

tr ans ( k , t- 1 ) =out ; 
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C(k,:) = [ y (2 .-length (y) ) y(l)]; 



end 

end 

% Normalize each row of C to 1 and 
% clip the data to 10^ (-6). 

psnorm=normod ( C , rows ) ; 
psnorm=clip (psnorm, 10^ (“6) ) ; 

% calculate the SPL of psnorm and store in M 
M=logl0 (psnorm) .*10; 

b. Kalman Filter using a Dynamic 0 

KALMOD.M 

This program filters the array of desired periodograms 
of the time series data using a kalman filter with a 
dynamic Beta. The entering arguments are; 

- <file name> of PSD array 

- sample size 

- sigma (already calculated if data was generated 
using "KFILTER.M" program. 

- total number of (n«m) points in the array (already 
calculated if "KFILTER.M" was used. 

- the kalman filter parameter "Beta" 

The subroutine outputs are; 

M= periodogram of filter output in dB 
C= periodogram of filter output in power 

% Define the function "KALMOD" 

function [M, C] =kalmod (psd , sampsize , sigma , totpts , betah) 

% KALMAN FILTER 

x=zeros (2 , sampsize) ; 
pointer=zeros( 2, sampsize) ; 

sv2=sigma^2 ; 

beta=betah 

betal=l 

% create loop for the number of rows in psd 
rows=totpts /sampsize ; 
for k=l:rows 

y=psd(k, : ) ; 
el=0. 0; 
e2=0.0; 
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dl=0 . 0 ; 



d2=0.0; 



x(l,l)=y(l) ; 
x(2,l)=y(l) ; 

tau=l. 0; 

% MAIN LOOP 

for t=l: (saropsize-l) ; 

kl=l. 0/ (tau+1. 0) ; 
xll=x(l,t)+kl*(y(t)-x(l,t) ) ; 
dll=dl-beta; 

ell=el+(l. 0/ (2.0*sv2) ) * ( (y (t+1) -xll) ^^2) ; 
cll=ell+dll; 

k2=l • 0 • 

xl2=x(l,t)+k2*(y(t)-x(l,t) ) ; 
dl2=dl+beta; 

el2=el+(l. 0/ (2. 0*sv2) ) * ( (y (t+1) -xl2) "2) ; 
Cl2=el2+dl2; 

kl=0 . 5 ; 

x21=x (2 , t) +kl* (y (t) -x(2 , t) ) ; 
d21=d2“beta; 

e21 = e2+(1.0/ (2.0*sv2) ) * ( (y (t+1) -x21) ''2) 
c21=e21+d21; 

k2=1.0; 

x22=x(2,t)+k2*(y(t)-x(2,t) ) ; 
d22=d2+beta; 

e22=e2+(1.0/(2.0*sv2) )*( (y(t+l)-x22) ^2) ; 
c22=e22+d22 ; 

% UPDATE STATES IN DYNAMIC PROGRAM. 

if clKc21 

x(l,t+l)=xll; 
el=ell ; 
dl=dll ; 
cl=el+dl ; 

pointer (1, t+1) =1 ; 

tau = tau+1; 
beta=betah ; 

else 

x(l,t+l)=x21; 
el=e21; 
dl=d21; 
cl=el+dl; 
tau=2 ; 

pointer ( 1 , t+l ) =2 ; 



90 



end 



beta=betal ; 



% 

% 



end 



if c22<cl2 

X (2 , t+1) =x22 ; 
e2=e22 ; 
d2=d22; 
c2=e2+d2 ; 
pointer (2 ,t+l) =2 ; 
beta=betal ; 

else 

x(2,t+l)=xl2; 
e2=el2 ; 
d2=dl2 ; 
c2=e2+d2 ; 
pointer ( 2 , t+1 ) =1 ; 
beta=betah ; 

end 



END MAIN LOOP 



BACKWARDS SMOOTHING AND SUBSTITUTION 
tau=l . 0 ; 



if cl<c2 

out=l ; 

else 



end 



out=2 ; 



y (sampsize) =x (out, sampsize) ; 

for t=sampsize:-l:2 

out=pointer (out, t) ; 
xout=x(out,t-l) ; 

if out==2 

tau=l . 0 ; 
y (t-1) =xout ; 

else 

tau=tau+l; 
y (t-1) =xout ; 

end 



trans (k, t-1) =out ; 

C(k,:)=[ y (2: length (y) ) y(l)]; 

end 

end 
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Normalize each row of C to 1 and 
clip the data to 10^ (-6). 

psnorm=normod ( C , rows ) ; 
psnorm=clip (psnorm, 10'' (-6) ) ; 

calculate the SPL of psnorm and store in M 
M=logl0 (psnorm) .*10; 

c. "Modified" Kalman Filter using Constant /5 
KALMANl.M 

This program filters the array of desired periodograms 
of the time series data using a ''modified” kalman filter 
with a constant Beta. The entering arguments are; 

- <file name> of PSD array 

- sample size 

- sigma (already calculated if data was generated 
using "KFILTER.M” program) . 

- total number of (n»m) points in the array (already 
calculated if "KFILTER.M” was used) . 

- the kalman filter parameter "Beta” 

The subroutine outputs are; 

M= periodogram of filter output in dB 
C= periodogram of filter output in power 

% Define the function "KALMANl” 

function [M, C]=kalmanl (psd, sampsize, sigma, totpts, beta) 

% KALMAN FILTER 

x=zeros (2 , sampsize) ; 
pointer=zeros ( 2 , sampsize) ; 

sv2=sigma''2 ; 

% create loop for the number of rows in psd 
r ows= t ot pt s / s amps i z e ; 
for k=l;rows 

y=psd(k, ; ) ; 

el=0.0; 

e2=0.0; 

dl=0.0; 

d2=0.0; 

x(l,l)=y(l) ; 
x(2,l)=y(l) ; 
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tau=l. 0; 



% MAIN LOOP 

for t=l: (sampsize-l) ; 

kl=l. 0/ (tau+1. 0) ; 
xll=x(l,t)+kl* (y (t) -x(l, t) ) ; 
dll=dl-beta ; 

ell=el+(1.0/ (2.0*sv2) ) * ( (y (t+1) -xll) ^2) ; 
cll=ell+dll; 

k2=l. 0; 

xl2=x(l,t)+k2*(y(t)-x(l,t) ) ; 
dl2=dl+beta; 

el2=el+(1.0/ (2.0*sv2) ) * ( (y(t+l) -xl2) '' 2 ) ; 
Cl2=el2+dl2; 

kl=0.5; 

x2l=x(2,t)+kl* (y (t) -x(2,t) ) ; 
d21=d2-beta; 

e21 = e2+(1.0/ (2.0*sv2) ) *( (y(t+l)-x21) "2) ; 
c21=e21+d21; 

k2=1.0; 

x22=x (2 , t) +k2* (y (t) -X (2 , t) ) ; 
d22=d2+beta; 

e22=e2+(1.0/ (2.0*sv2) ) * ( (y (t+1) -x22) ^2) ; 
c22=e22+d22 ; 

% UPDATE STATES IN DYNAMIC PROGRAM. 

if C1KC21 

x(l, t+1) =xll; 
el=ell ; 
dl=dll ; 
cl=el+dl ; 

pointer (1, t+1) =1 ; 

tau = tau+1; 

else 

x(l,t+l)=y (t) ; 
el=e21; 
dl=d21; 
cl=el+dl; 
tau=2 ; 

pointer ( 1 , t+1 ) =2 ; 

end 

if c22<cl2 

x(2,t+l)=y(t) ; 
e2=e22 ; 
d2=d22; 
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% 

% 



end 



c2=e2+d2 ; 
pointer ( 2 , t+1 ) =2 ; 

else 

x(2,t+l)=xl2; 
e2=el2 ; 
d2=dl2 ; 
c2=e2+d2 ; 
pointer ( 2 , t+1 ) =1 ; 

end 



END MAIN LOOP 

BACKWARDS SMOOTHING AND SUBSTITUTION 



tau=l . 0 ; 
if cKc2 



else 

end 



out=l ; 
out=2 ; 



y (sampsi 2 e)=x(out, sampsize) ; 

for t=sampsize:-l:2 

out=pointer (out,t) ; 
xout=x(out, t-1) ; 

if out==2 

tau=l . 0 ; 
y (t-l)=xout; 

else 

tau=tau+l ; 
y (t-l)=xout; 

end 



trans (k , t-1 ) =out ; 

C(k,:)=[ y (2: length (y)) y(l)]; 

end 

end 

% Normalize each row of C to 1 and 
% clip the data to 10'' (-6). 

psnorm=normod(C,rows) ; 
psnorm=clip(psnorm, 10'' (-6) ) ; 

% calculate the SPL of psnorm and store in M 
M=logl0 (psnorm) . *10 ; 
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d. Kalman Filter using a Constant /3 



KALMAN. M 

This program filter the array of desired periodograms 
of the time series data using a kalman filter with a 
constant Beta. The entering arguments are; 

- <file name> of PSD array 

- sample size 

- sigma (already calculated if data was generated 
using ''KFILTER.M” program) . 

- total number of (n*m) points in the array (already 
calculated if "KFILTER.M" was used) . 

- the kalman filter parameter "Beta" 

The subroutine outputs are; 

M= periodogram of filter output in dB 
C= periodogram of filter output in power 

% Define the function "KALMAN" 

function [M, C]=kalman (psd, samps ize, sigma , totpts, beta) 

% KALMAN FILTER 

x=zeros (2 , sampsize) ; 
pointer=zeros (2 , sampsize) ; 

sv2=sigma^2 ; 

% create loop for the number of rows in psd 
rows=totpts/sampsize ; 
for k=l;rows 

y=psd(k, : ) ; 
el=0. 0 ; 
e2=0. 0; 

dl=0 . 0 ; 
d2=0.0; 

x(l,l)=y(l) ; 
x(2,l)=y(l) ; 
tau=l. 0 ; 

% MAIN LOOP 

for t=l: (sampsize-1) ; 

kl=1.0/ (tau+1.0) ; 
xll=x(l,t)+kl*(y(t) -x(l,t) ) ; 
dll=dl“beta ; 

ell=el+(1.0/ (2.0*sv2) ) * ( (y (t+1) -xll) ^2) ; 
cll=ell+dll; 
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k2=1.0; 

xl2=x(l, t) +k2* (y (t) “X(l, t) ) ; 
dl2=dl+beta; 

el2=el+(1.0/(2.0*sv2) )*( (y(t+l)-xl2) "2) ; 
cl2=el2+dl2 ; 

kl=0 . 5 ; 

x21=x(2 ,t) +kl* (y (t) -x(2 ,t) ) ; 
d21=d2-beta; 

e21 = e2+(1.0/(2.0*sv2))*((y(t+l)-x21)^2) 
c21=e21+d21; 

k2=1.0; 

x22=x(2,t) +k2* (y(t)-x(2,t)) ; 
d22=d2+beta ; 

e22=e2+ (1. 0/ (2 . 0*sv2) ) * ( (y (t+1) -x22 ) ^2) ; 
C22=e22+d22; 

% UPDATE STATES IN DYNAMIC PROGRAM. 

if clKc21 

x(l,t+l)=xll; 
el=ell ; 
dl=dll ; 
cl=el+dl ; 

pointer (1, t+1) =1 ; 

tau = tau+1; 

else 

x(l,t+l)=x21; 
el=e2l; 
dl=d2l; 
cl=el+dl; 
tau=2 ; 

pointer ( 1 , t+l ) =2 ; 

end 

if c22<cl2 

x(2, t+l)=x22 ; 
e2=e22 ; 
d2=d22 ; 
c2=e2+d2 ; 
pointer ( 2 , t+l ) =2 ; 

else 

x(2,t+l)=xl2; 
e2=el2 ; 
d2=dl2 ; 
c2=e2+d2 ; 
pointer ( 2 , t+l ) =l ; 
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% 



END MAIN LOOP 



% BACKWARDS SMOOTHING AND SUBSTITUTION 



tau=l . 0 ; 



if cl<c2 

out=l ; 

else 



end 



out=2 ; 



y(sainpsize)=x(out,sainpsize) ; 

for t=sainpsize:-l:2 

out=pointer (out, t) ; 
xout=x (out, t-1) ; 

if out==2 

tau=l . 0 ; 
y (t-1) =xout ; 

else 

tau=tau+l; 
y (t-1) =xout ; 

end 

trans (k , t-1 ) =out ; 

C(k,:)=[ y(2:length(y) ) y(l)]; 

end 

end 

% Normalize each row of C to 1 and 
% clip the data to 10^ (-6). 

psnorm=normod ( C , rows ) ; 
psnorm=clip(psnorm, 10^ (-6) ) ; 

% calculate the SPL of psnorm and store in M 
M=logl0 (psnorm) . *10 ; 

2 . Supporting Matlab Code 

a. Signal Generator 

% KFILTER.M R.C. Adamo 

% The purpose of this program is to connect all ap- 
% propriate subroutines allowing the user to generate 
% a single stationary frequency component test signal 
% for processing. The user is prompted for; 

% 1) Input frequency 
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% 2) Sampling frequency 

% 3) Sample size (Power of 2 please) 

% 4) Signal length (Power of 2 please) 

% 5) Input SNR 

% The signal is processed in non-overlapping 
% segments zero-padded to twice the sample size. 

% The following subroutines are called; 

% 1) PER.M 

% 2) KALMODl.M 

% (KALMAN. M, KALMANl.M AND KALMOD.M can also 

% be used if desired. User must delete the 

% '•%" preceding the desired subroutine call-up 

% in this program) 

% After processing, the square of the magnitudes of the 
% FFT's are saved in arrays of (signal length/ sample 
% size) rows by sample size columns as; 

% - rpsd.dat 

% - hwpsd.dat 

% - kpsd_2B_mod.dat (if you modify this program to 

% run the other kalman algorithms you must save 

% that data manually) 

% For a signal length of (128*128 ) , allow 10 mins per Kalman 
% filter utilized for total run time. 

clear 

% Input desired frequency 

fl=input( 'Enter Signal freq less than 256 Hz ') 

% Input sampling frequency 

fs=input( 'Enter sampling freq ') 

% Input sample size 

size=input( 'Enter sample size (power of 2 pis) ') 

% Input the number of points to be generated. 

totpts=input (' Enter total # of pts (power of 2 pis)') 
no=0:l: (totpts-1) ; 

% Input the desired SNR of this signal 
snr=input ( 'Enter SNR in dB ') 
sigma=sqrt( (10^ (-snr/10) ) *. 5) 

% Generate sinusoid only 
thetal=2*pi*f l*no; 
yl=sin (thetal/fs) ; 
ycomp=yl ; 

% Add Gaussian noise 
rand( 'normal' ) 
rand ( ' seed' , 5) 

ydirt=ycomp+sigma*rand(ycomp) ; 
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% define beta parameters 

bl=input( ' Enter desired filtering beta ') 

% Call subroutines to compute rec,rpsd,h,hwpsd,kw,kwpsd 
[rec,rpsd,h,hwpsd,freq]=per(ydirt,size, fs) ; 

% [k_lB,kpsd_lB]=kalman(rpsd, size, sigma, totpts,bl) ; 

% [k_lB_mod,kpsd_lB_mod]=kalmanl(rpsd, size, sigma, totpts,bl) ; 
% [k_2B,kpsd_2B]=kalmod (rpsd, size, sigma, totpts , bl) ; 

[k_2B_mod,kpsd_2B_mod]=kalmodl (rpsd, size, sigma , totpts , bl) ; 

save 

save rpsd.dat rpsd /ascii 

save hwpsd.dat hwpsd /ascii 

save kpsd_2B_mod.dat kpsd_2B_mod /ascii 

quit 



b. Subroutine "PER.H" 



% PER.M R.C. Adamo 

% Compute the rectangular window and the hamming window 
% periodograms and the frequency vector for the desired 
% signal. The input arguments are signal name, the 

% sample size of the transforms and the sampling frequency. 

% The data segments are zero-padded to twice the sample 
% size. 

% The output data of the subroutine are: 

% M=periodogram of unwindowed time series in dB 

% B=periodogram of unwindowed time series in power 

% C=periodogram of hamming windowed time series 

% in dB 

% D=periodogram of hamming windowed time series 

% in power 

% E=sampling frequency vector used for plotting 

% individual periodograms. 

% Define the function "PER” 

function [M,B, C,D, E]=per (signal, sampsize, sampfreq) 

rows=length (signal) /sampsize; 
for k=l trows 

% zero pad sample size pt data sets to twice the sample 
% size. 

temp (1: sampsize) =signal( (k-1) *sampsize+l ;k*sampsize) ; 
temp (sampsize+1; 2*sampsize)=zeros (1: sampsize) ; 

% Square the magnitude of the fft of the padded sequence 
ztemp=(abs (f ft (temp) ) ) . ^2 ; 

% Store the left half of ztemp as Matrix where 
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% Matrix is n by sample size matrix 
B(k, :)=ztemp(l:sampsize) ; 

clear temp ztemp 
end 

% Normalize each row of B and call it psnorm 
psnorm=normod(B,rows) ; 

% Clip psnorm to minimum value of 10''-6 
psnorm=c lip (psnorm, (-6) ) ; 



% calculate the SPL of psnorm and store in M 
M=loglO (psnorm) .*10; 

[C,D]=hamm( signal, sampsize) ; 

E=fvec(sampfreq, sampsize) ; 

c. Hamming Windotf Subroutine 

% HAMM.M R.C. Adamo 

% Compute the periodogram of desired signal after applying 
% a Hamming window to each time series data sample of 
% defined sample size. The input arguments are the signal 
% name and the transform sample size. The data segments are 
% zero-padded to twice the sample size. 

% The output data of the subroutine are: 

% M=periodogram of hamming windowed time series 

% in dB 

% B=periodogram of hamming windowed time series 

% in power 

% Define the function "HAMM" 

function [M, B]=hamm(signal, sampsize) 

% Create hamming window of in-putted sample size 
w=hamming (sampsize) ; 

rows=length (signal) /sampsize; 
for k=l crows 

% apply hamming window and zero pad data to twice the sample 
% size. 

temp (1: sampsize) =w' .*signal( (k-1) *sampsize+l:k*scunpsize) ; 
temp ( sampsize+1 : 2*sampsize) =zeros ( 1 : sampsize) ; 

% Square the fft magnitude of the padded sequence 
ztemp=(abs(fft(temp) ) ) . ^2; 
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% store the left half of ztemp as Matrix where Matrix 
% is n by sample size; n=totpts/ sample size. 

Matrix (k, : )=ztemp(l:sampsize) ; 

clear temp ztemp 
end 

% Normalize each row of Matrix and call it psnorm 
psnorm=normod (Matrix, rows) ; 

B=Matrix; 

% Clip psnorm to minimum value of 10^-6 
psnorm=clip(psnorm, 10^ (”6) ) ; 

% Calculate the SPL of psnorm and store in M 
M=logl0 (psnorm) .*10; 

d. Subroutine "NORMOD.M" 

% NORMOD.M R. C. Adamo 

% This program normalizes a matrix such 

% that the largest value in each row is 

% 1 and the other values are scaled relative to 1. 

% The entering arguments are the matrix name and 

% the number of rows in the matrix. 

% Define the function "NORMOD" 
function M=normod (matrix, rows) 

for k=l:rows 

M(k, : )=(l/max(matrix(k, : ) ) ) * matrix(k,:); 

end 



e. Subroutine "CLIP.M" 

% CLIP.M R.C. ADAMO 

% This subroutine clips any matrix to some 
% minimum value. The entering arguments are the 

% matrix to be clipped and the minimum value the user 

% wishes to clip it to. 

% Define the function "CLIP" 

function M=clip (matrix, clip) 

index = matrix < (clip) ; 
corr = index * clip; 
corfac = 1 - index; 
adjust = corfac. *matrix; 

M = corr + adjust; 
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f. Subroutine " NORMAL. M" 



% NORMAL. M R. C. Adamo 

% This program normalizes a matrix such 
% that the largest value in the matrix are 
% 1 and the other values are scaled relative to 1. 

% Define the function "NORMAL” 
function M=normal (matrix) 

M= ( 1/max (max (matrix) ) ) *matrix; 

g. Subroutine "FVEC.M" 



% FVEC.M R.C. Adamo 

% CREATE THE FREQUENCY VECTOR USED IN PLOTTING 
% A PERIODOGRAM. fs IS THE SAMPLING FREQUENCY 
% AND X IS THE SAMPLE SIZE. 

function f=fvec(fs,x) 

n=2*x; 

f=fs* (0:n-l) /n; 

h. Grey-tone imaging FORTRAN 77 Code 

C LOFAR.F R.C. ADAMO 

C This program reads and displays processed signals 
C in a LOFAR presentation. It is used to compare 
C the performances of different processing techniques. 

C The user must edit this program to ensure input filedef 
C is properly defined. Also this program must be edited to 
C add a proper title on output. 

C 

C OPEN GKS, OPEN WORKSTATION OF TYPE 1, ACTIVATE WORKSTATION 
C 

CALL GOPKS (6,IDUM) 

CALL GOPWK (1, 2, 1) 

CALL GACWK (1) 

C 

C INVOKE DEMO DRIVER 
C 

call thafto(ierr) 

C 

C DEACTIVATE AND CLOSE WORKSTATION, CLOSE GKS. 

C 

CALL GDAWK (1) 

CALL GCLWK (1) 

CALL GCLKS 
STOP 
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END 



SUBROUTINE THAFTO ( TERROR) 

C 

C USAGE 
C 

C ARGUMENTS 
C ON OUTPUT 
C 
c 
c 
c 

C I/O 

c 
c 
c 
c 
c 
c 
c 

C PRECISION 
C 

C LANGUAGE 
C 

C REQUIRED ROUTINES HAFTON 

C 

C REQUIRED GKS LEVEL OA 

C 

C Z contains the values to be plotted. 

C 

REAL Z(128,128) 

INTEGER ROWS, SIZE 

SAVE 

Specify coordinates for plot titles. The values TX and TY 
define the center of the title string in a 0. to 1. range. 

DATA TX/0.15/, TY/0.9469/ 

Specify low (FLO) and high (FHI) contour values, and NLEV 
unique contour levels. NOPT determines how Z maps onto the 
intensities, and the directness of the mapping. 

DATA FLO/0.0/, FHI/0.0/, NLEV/1/, NOPT/1/ 

Initialize the error indicator 

TERROR = 0 

Define files to be read in. 
open(l, file=' rpsd.dat' ) 



CALL THAFTO (TERROR) 



TERROR 

An integer variable 

=0, if the test was successful, 

= 1, the test was not successful. 

If the test is successful, the message 
HAFTON EXECUTED — SEE PLOTS TO CERTIFY 
is printed on unit 6. In addition, 1 half-tone 
frame is produced on the machine graphics 
device. In order to determine if the test 
was successful, it is necessary to examine 
the plot. 

Single 
FORTRAN 77 
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C open(l,file='h_snrcompl.dat' ) 

C open ( 1 , f i le= ' k_2B_mod_13 5k . dat ' ) 

C 

C Input data array dimensions 

print *, ' Enter array size as rows, column ' 
read *, rows, size 
rows=128 
size=128 

Read in the desired data to be plotted 
DO 20 1=1, rows 

READ(1,*,END=999) (Z(J,I) ,J=l,size) 

20 CONTINUE 

Select normalization trans 0 for plotting title. 

CALL GSELNT (0) 

Call WTSTR to write the plot title. 

Plot upper title for output 
CALL WTSTR (TX,TY, 

1 'LOFAR GRAM; RECT WINDOW; Base SNR out=9 

1 dB',2,0,-1) 

1 'LOFAR GRAM; HAMM WINDOW; Base SNR out=9 

1 dB',2,0,-1) 

1 'LOFAR GRAM; ''Mod” Kalman dynamic Beta=135k' 

1 , 2 , 0 ,- 1 ) 

Print frequency scale on output 
CALL WTSTR (0.10,0.08, 

1 '0 10 20 30 
1 ',2, 0,-1) 

CALL WTSTR (0.40,0.055, 

1 'frequency Hz ',2, 0,-1) 

Entry HAFTON allows user specification of plot parameters. 

CALL HAFTON ( Z , size , size, rows , FLO , FHI , NLEV, NOPT ,0,0,0.) 
CALL FRAME 
C 

WRITE (6,1001) 

RETURN 

1001 FORMAT (' HAFTON TEST EXECUTED — SEE PLOTS TO CERTIFY') 
C 

999 write(6,llll) 

1111 format ( 'unexpected end of file on input: program 
C 1 terminated') 

END 
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i. Subroutine "NAMES .M" 



% NAMES. M R.C. Adamo 

% Create the titles used for the comparison plots 

dBname=[ sprint f ('SNR=%g' ,snr) , 'dB, ' ] 

Betaname=sprintf ( ' BETA=%g ' , bl ) 
freqname=[ sprint f ( ' freq=%g' , f 1) , 'Hz' ] 

Rectitle=[ 'REC Win, ' ,dBname, freqname] 

HWtitle=[ 'HAM Win, ' ,dBname, freqname] 

Kaltitle= [ 'Mod Kal, dynamic beta, ',Betaname] 

Kltitle =['Mod Kal, constant beta, ',Betaname] 

K2title =['Kal Out, dynamic beta, ',Betaname] 

K3title =['Kal Out, constant beta ,',Betaname] 
Betahigh=sprintf ( 'Beta high = %g',bl) 

Betalow= sprintf ( ' Beta low = %g',l) 

j . Plotting Subroutine 

% PLOTT.M R.C. Adamo 

% This subroutine plots and compares the periodograms 
% of the rectangular and Hamming windows and the Kalman 
% filter outputs in dB. If other than the default 
% kalman algorithm is used, the user must modify line 18 
% of this program to plot the desired kalman output. 

% This program must be run after executing the program 
% "names” in matlab. 

for i=l ; 1 : totpts/size 

clg 

subplot (221) 

plot (freq(l: size) ,rec(i, : ) ) 

title (Rectitle) 

xlabel (' frequency Hz') 

ylabel ( 'magnitude dB') 

plot (freq(l: size) ,h(i, ; ) ) 

title (HWtitle) 

xlabel (' frequency Hz') 

ylabel ( 'magnitude dB') 

plot (freq(l: size) ,k_2B_mod(i, : ) ) 

title (Kaltitle) 

xlabel (' frequency Hz') 

ylabel ( ' magnitude dB ' ) 

text ( . 60 , . 30, Betahigh, ' sc' ) 

text ( . 60, . 25, Betalow, ' sc' ) 

text (. 60 ,. 20, sprintf ( 'Noise realization %g',i),'sc') 
pause 
end 
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APPENDIX B. PERFORMANCE RELATIVE TO SPECTRAL RESOLUTION 



The ability of the Kalman filter to detect two components 
at nearby spectral locations was examined. Three test cases 
were evaluated to determine the minimum spectral separation 
needed by the Kalman filter to distinguish between two 
components. The same data sets were also processed using the 
Hamming window for comparison. 

All three test cases consisted of three sinusoidal 
components in additive Gaussian noise. The frequencies in 
test case 1 were 10.7 Hz, 11.2 Hz and 20.1 Hz at input of 
-9 dB, -12 dB and -15 dB, respectively. This case provides a 
2 bin separation between fl (10.7 Hz) and f2 (11.2 Hz). This 
signal was processed using the Hamming window and the Kalman 
('•modified” dynamic /3=135k) filter. Test frequencies in test 
case 2 were 10.7 Hz, 11.45 Hz and 20.1 Hz at input of -9 
dB, -12 dB and -15 dB (3 bin separation between fl and f2) , 
respectively. Finally, test frequencies in test case 3 were 
10.7 Hz, 11.7 Hz and 20.1 Hz at input of -9 dB, -12 dB 
and -15 dB (4 bin separation between fl and f2) , respectively. 
Figure 42 thru Figure 44 shows the LOFAR outputs for the 
Hamming window and the Kalman ("modified” dynamic j8=135k) for 
cases 1, 2 and 3, respectively. The mainlobe broadening of 
the Hamming window over that of the Kalman filter is obvious 
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in all three cases. In addition, the Kalman filter clearly 
detects the distinct frequencies of fl and f2 in case 3 (4 bin 



separation) where the Hamming window process can not. 
difficult to pull out f3, the -15 dB component, 
either the Hamming or Kalman filter in any of the three 



It is 
using 
cases . 
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Figure 42. 2 BIN SEPARATION. Sinusoid (fl=10.7Hz at SNRi,= 

-9dB, f2=11.2Hz at SNRi=-12dB and f3=20.1Hz at SNR^=~15dB) 

comparing Hamming vs Kalman ("modified" dynamic /3=135k) . 
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Figure 43. 3 BIN SEPARATION. Sinusoid (fl=10.7Hz at SNR^=- 
9dB, f2=11.45Hz at SNR^=-12dB and f3=20.1Hz at SNRi„=-15dB) 
comparing Hamming vs Kalman ("modified” dynamic /S=135k) . 
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Figure 44 . 4 BIN SEPARATION. Sinusoid (fl=10.7Hz at 

9dB, f2=11.7Hz at SNRi=~12dB and f3=20.1Hz at SNi?i„=-15dB) 
comparing Hamming vs Kalman ("modified" dynamic /3=135k) . 
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APPENDIX C. EFFECTS OF DYNAMIC SPECTRAL COMPONENTS 

LOFAR outputs are extremely useful in providing the 
additional processing gain needed to detect signal in noise at 
SNR^ of 6 dB or less. Up to this point the experimental data 
was comprised of stationary frequency components. It has been 
demonstrated that the human eye is well suited to linearly 
regress the data displayed via the LOFAR output for stationary 
frequency components. The effects that dynamic spectral 
components will have on this capability is examined in this 
appendix. The time series data for this experiment was 
generated with three frequency components. One component is 
a stationary frequency of 20.45 Hz at an SNR^, of -15 dB. The 
second spectral component is a linear ramp frequency with an 
SNRj„ of -12 dB starting at 11.45 Hz and ending at 31.45 Hz. 
The third component is a dynamic component with an SNRj„ of -9 
dB starting and ending at 15.7 Hz and completing one sine wave 
cycle of amplitude of 2.5 Hz. 

Figure 4 5 is the LOFAR output of this data filtered by the 
Hamming window and the Kalman ("modified" dynamic |3=135k) 
filter. Both dynamic components are clearly visible in both 
the Hamming and Kalman LOFAR outputs. LOFAR output allows the 
humman eye to observe the structure of dynamic spectral 
components at SNR's^ down to low as 6 dB. 



Ill 





Figure 45. LOFAR comparison of Hamming window vs Kalman 
(••modified” DYN /3=135k) of dynamic freqs (sine wave at 
SNRi„=-9 dB, ramp at SNRj„=-12 dB, fixed at SNRi„=~15 dB) . 



112 




APPENDIX D. LOFAR OUTPUT USERS GUIDE 



A. USERS GUIDE 

This appendix provides the computer code and the step-by- 
step instructions necessary to produce a LOFAR output of any 
array of data. The user is required to load the data of 
interest into the same directory that the LOFAR code "image, f" 
is using. The data file must be named "input.dat". When the 
data is in place the following procedures should be followed: 

i) For output to be displayed on the Sun workstation; 

1. Rename data file to "input.dat". 

2. At the prompt type "ncargf77 -o image image.f". 
This compiles the program and names the executable 
file "image". 

3. At the prompt type "ncargrun -o image. cgm image". 
This executes image and saves the LOFAR output as 
" image . cgm" 

4. The program prompts the user to input the number of 
rows and number of columns in their data. At the 
prompt type "(# of rows) <space> (# of columns)". 
The program then reads in the data and builds the 
LOFAR output. 

5. At the prompt type "ctrans -d sunview image. cgm". 
This translates the "image. cgm" to the appropriate 
form for display on the Sun workstation. 

ii) Hard copy LOFAR output 

Once the procedures described above have been executed, 
hard copy of the LOFAR output ("image. cgm") can be produced on 
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the tektronix 4014 (tek4014) terminal in the computer lab Sp- 
3 01 if desired. To obtain a hard copy output follow these 
procedures : 

1. Using the tek4014 terminal in Sp-301 log on to the 
VAX and type ''tek4014” at the prompt <term(vt200) >. 

2 . Remote logon to the DSP used to produce the LOFAR 
output on. Again type "tek4014” at the 
<term (vt200) > prompt. 

3. Change to the directory where the LOFAR output file 
” image . cgm” resides. 

4. Type "ctrans -d t4010 image. cgm". 

5. It takes between 1-2 minutes to display the LOFAR 
output on the tektronix. Once plotting is complete 
a hard copy can be obtained by pressing the "copy 
button" on the Tektronix 4631 Hard Copy Unit. 



B. COMPUTER CODE 

C IMAGE. F R.C.ADAMO 

C This program reads and displays data sets of 
C "n by m" dimensions. The user must ensure the data 
C is saved under the name of "INPUT.DAT". The user is 
C prompted for the dimensions of the data set. The 
C number of rows and columns must be integers. 

C 

C OPEN GKS, OPEN WORKSTATION OF TYPE 1, ACTIVATE WORKSTATION 
C 

CALL GOPKS (6,IDUM) 

CALL GOPWK (1, 2 , 1) 

CALL GACWK (1) 

INVOKE DEMO DRIVER 

call thafto(ierr) 

DEACTIVATE AND CLOSE WORKSTATION, CLOSE GKS. 

CALL GDAWK (1) 

CALL GCLWK (1) 

CALL GCLKS 
STOP 
END 
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SUBROUTINE THAFTO (TERROR) 



C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 



USAGE 



CALL THAFTO (lERROR) 



ARGUMENTS 

ON OUTPUT I ERROR 

An integer variable 

=0, if the test was successful, 

= 1, the test was not successful. 

I/O If the test is successful, the message 

HAFTON EXECUTED — SEE PLOTS TO CERTIFY 

is printed on unit 6. In addition, 1 half-tone 
frame is produced on the machine graphics 
device. In order to determine if the test 
was successful, it is necessary to examine 
the plot. 

PRECISION Single 

LANGUAGE FORTRAN 77 



REQUIRED ROUTINES HAFTON 

REQUIRED GKS LEVEL OA 



Z contains the values to be plotted. 

REAL Z(128,128) 
INTEGER ROWS, SIZE 
SAVE 



Specify coordinates for plot titles. The values TX and TY 
define the center of the title string in a 0. to 1. range. 

DATA TX/0.15/, TY/0.9469/ 

Specify low (FLO) and high (FHI) contour values, and NLEV 
unique contour levels. NOPT determines how Z maps onto the 
intensities, and the directness of the mapping. Execute 
"man HAFTON” for a detailed description of these variables. 

DATA FLO/0.0/, FHI/0.0/, NLEV/1/, NOPT/1/ 

Initialize the error indicator 



TERROR = 0 
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Define file to be read in. 
open (1, f ile=' : iput.dat') 



Is data file called INPUT.DAT? 



print 
print 
print 
print 
print 
print 
read *, k 

if (k.eq.2) go to 991 



Is your data file named INPUT.DAT?' 
Enter 1 for yes and ' 

2 for no. ' 

/ 

If you answer no the program will ' 
terminate. Rename file and try again. ' 



Input data array demens ions 

print *, ' Enter array size as rows, column ' 
read *, rows, size 



Read in the desired data to be plotted 



DO 20 1=1, rows 

READ(1,*,END=999) (Z(J,I) , J=l,size) 

20 CONTINUE 



Select normalization trans 0 for plotting title. 
CALL GSELNT (0) 



Call WTSTR to write the plot title. 

CALL WTSTR (TX,TY, 

1 'Your Plot tilte goes here' , 2 , 0, -1) 

CALL WTSTR (0.15,0.08, 

1 'Scale data could go here' ,2, 0,-1) 



Entry HAFTON allows user specification of plot parameters. 

CALL HAFTON ( Z , size , size, rows , FLO, FHI , NLEV, NOPT, 0,0,0.) 
CALL FRAME 
C 

WRITE (6,1001) 

RETURN 

C 

1001 FORMAT (' HAFTON TEST EXECUTED — SEE PLOTS TO CERTIFY') 
C 

999 write(6, 1111) 

1111 format ( 'unexpected end of file on input: program 
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1 terminated') 
go to 9999 

991 print *, 'rename input file to INPUT.DAT and try again' 
9999 END 
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